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Introduction 

In  this  report  we  describe  work  completed  under  Subcontract  No. 
N68305-75-C-004  (Proposal  UCB-Eng.-3889).  As  we  frequently  refer  to  our 
previous  work  [1-5],  it  will  prove  helpful  to  the  reader  to  be  somwhat 
familiar  with  these  references. 

In  Section  I we  present  theory  and  algorithms  for  large  displacement 
contact-vimpact  analysis  in  two  dimensions(i .e.  plane  stress,  plane  strain 
and  axi symmetric) . This  work  builds  upon  earlier  developments  documented  in 
[2]  and  [3]  and  represents  the  completion  of  our  theoretical  work  in  this 
area.  The  theory  encompasses  a wide  range  of  contact- impact  problems  and 
allows  for  a completely  arbitrary  contact  surface  development,  stick,  slip 
and  frictional  sliding  conditions,  and  impact-release  conditions  covering 
the  full  range  of  contact  possibilities. 

In  Section  II  we  describe  some  new  developments  regarding  the  Hertzian 
algorithm,  which  has  been  extensively  documented  in  previous  publications; 
see  [1-5], 

Section  III  contains  sample  problems  which  employ  the  algorithms 
described  in  Section  I and  also  some  studies  involving  the  Hertzian  algorithm. 

Anticipating  the  need  for  an  efficient  shell  element  for  crash 


configuration  modelling,  we  have  performed  a pilot  study  of  a plate  bending 
element.  The  results  of  this  study,  which  are  encouraging,  are  contained 
in  Section  IV. 
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I.  Large  Displacement  Contact-Impact  Theory 
1 . Introduction 

In  the  following  section^  we  present  a large  displacement  contact-impact 
theory.  In  Section  1-2  we  set  notations  and  establish  the  structure  of  the 
local  equations  for  a typical  contactor  node  and  target  element  boundary. 

The  algorithm  for  handling  stick,  slip  and  frictional  contact  conditions  is 
discussed  in  Section  1-3.  This  strategy  is  sufficient  for  the  quasi-static 
problem  and  governs  the  iterations  within  a time  step  in  a dynamic  problem. 
The  updating  ir  dynamic  problems  manifested  by  discrete  impact-release 
conditions  is  described  in  Section  1-4. 
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2.  Local  Structure  of  the  Equations 

To  explicate  the  structure  of  the  equations  we  wish  to  solve,  we  shall 
consider  in  this  section  an  idealization  consisting  of  one  contactor  node 
and  one  target  element  boundary  (see  Fig.  1-1).  Quantities  associated  with 
the  contactor  node  will  possess  a subscript  k,  and  those  associated  with  the 
target  nodes  will  possess  a subscript  £ or  £+1.  Throughout  this  section  we 
shall  concern  ourselves  only  with  equations  which  pertain  to  these  three 
nodes.  In  two-dimensional  analysis  (i.e.  plane  stress,  plane  strain  or 
axi symmetric)  there  are  two  displacement  degrees-of- freedom  associated  with 
each  of  the  three  nodes  and,  in  addition,  two  contact  force  degrees -of- 
freedom  associated  to  the  contactor  node;  a total  of  eight  degrees-of- 
freedcm  for  the  three  nodes  considered  here.  Thus  there  are  six  equations 
of  motion  to  be  satisfied  and,  if  the  bodies  are  in  contact,  two  conditions 
of  compatibility.  Consider  the  case  in  which  the  contactor  node  is  in 
contact  with  the  target  element  boundary.  Then  if  a is  the  nondimensional 
location  parameter  for  the  target  element  boundary,  defined  by  (see  Fig.  1-2): 

° = ^xlk  " xl^2  + (x2k  " x2/31/2  /L  * 

L = ^xi,  £+1  ‘ xn)2  + (x2,  £+i  ' x2£>231/2  ’ 

the  eight  equations  to  be  satisfied  for  the  three  nodes  are: 
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the  compatibility  conditions.  The  superscript  indicates  to  which  body  the 

quantity  pertains  and  the  subscript  y indicates  the  coordinate  direction 

z^,  y = 1,2.  Note  that  the  components  K and  K are  functions  of  the  dis- 
1 2 

placements  of  B and  B , respectively. 

The  compatibility  conditions  amount  to  the  second  and  sixth  equations 

in  1-1.  When  the  bodies  are  not  in  contact  these  conditions  are  ignored 

and  x . = o. 
yk 

The  solution  of  the  matrix  equations  is  developed  by  employing  a 
temporal  discretization  which  results  in  a nonlinear  algebraic  problem  to  be 
solved  at  each  time  step.  The  Newmark  family  of  algorithms  is  employed 
by  us  to  temporally  discretize  the  equations  and  we  confine  our  attention 
to  implicit  methods.  This  aspect  of  our  work  has  oeen  described  previously 
(see  e.g.  [5])  and  we  shall  not  repeat  the  details  here.  The  resulting 
temporally  discretized  system,  for  the  nodes  in  question,  becomes: 
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where  the  subscripts  n and  n+1  indicate  the  time  step  at  which  the  quantity 
is  evaluated.  At  - tn+1  - tR  is  the  length  of  the  time  interval,  and  8 is 
the  Newmark  parameter. 

To  solve  the  nonlinear  algebraic  problem  at  each  time  step,  the  Newton- 
Raphson  method  is  employed.  The  linear  equations  used  in  this  procedure, 
for  the  nodes  in  question,  are  given  as  follows: 
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where  the  superscript  in  parentheses  indicates  the  iteration  number  and 
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denotes  the  tangent  stiffness  matrix.  The  formation  of  the  tangent  stiffness 
follows  the  standard  rules  {in  particular,  in  the  case  of  linear  elastic 
bodies  it  is  the  usual  stiffness  matrix).  The  iterative  process  is  defined 
by 

w(°)  _ w 

-n+1  " -n  ’ 


w(i+l)  = W(D  + Aw(i) 

-n+1  -n+1  %+l 


and  when  w^j  satisfies  a convergence  test,  for  some  i,  then  wn+-|  ^ 
Here  w is  defined  by 
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is  called  the  contact  stiffness.  Specification  of  it. 
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L*2. 


n+1 


are  the  unique  aspects  of  the  contact-impact  algorithm.  In  the  next  section 
we  describe  the  way  this  is  done. 

In  passing,  we  note  that  formally  setting  m to  o provides  an  algorithm 
for  the  quasi-static  case. 
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3.  Computing  Strategy  Within  a Time  (Load)  Step  * 

In  this  section  we  describe  the  strategy  involved  in  solving  the  non-  jj 

i 

* 

linear  algebraic  contact  problem  within  a time  step  (or,  equivalently,  within 
a load  step  in  the  quasi-static  case).  The  unknowns  involved,  as  described 
in  the  previous  section,  are  the  displacements  and  contact  forces  at  the  end 
of  the  step.  Treatment  of  velocities  and  accelerations,  and  updating  of 
contact  forces,  when  impact  and  release  effects  are  present,  are  considered 
in  the  following  section. 

For  consideration  of  sliding  contact,  it  is  convenient  to  work  in  a 
coordinate  system  naturally  defined  by  the  target  segment.  A system  of  this  ’ 

kind  can  be  constructed  by  aligning  the  coordinate  axes  in  the  tangential 

i 

(s)  and  normal  (n)  directions  to  the  target  element  boundary,  with  origin 
located  at  node  £ (see  Fig.  1-3).  If  e denotes  the  angle  between  the  s and 

l 

Z-.  axes,  measured  counterclockwise  from  z-,,  then  \ 

i 

8 = arc  tan  [U2>  t+,  - £+,  - x,t)]  . j 

and  vectors  nay  be  resolved  in  the  usual  way  into  tangential  and  normal 
components,  e.g. 
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where  c = cos  e and  s = sin  6.  Thus  the  vectors  in  the  s,  n - system 
corresponding  to  and  ^ are 

is  " (,sk  0 -(,-“>'sx  ■“  Tsk)T  and 

in  ' (tnk  0 -(1-“,Tnk  Tnk)T  • 


respectively,  and  likewise  for  x and  x . 
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At  the  end  of  each  iteration,  one  of  three  conditions  can  hold  for  a 
typical  contactor  node  k;  the  condition  is  determined  by  the  contactor  code 
V ^ ^k  = °*  then  contactor  node  k not  1n  contact;  if  ik  > o,  then 
node  k is  in  contact.  The  code  ik  = 1 signifies  the  stick  condition  and 
ik  = 2 signifies  the  sliding  condition.  The  contact  code  determines  what 
is  assembled  into  the  contact  stiffness  and  the  right-hand  side  vectors  t 
and  x.  Specifically,  we  have  the  following  situations: 


ik  = °: 
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and  fd  is  the  dynamic  coefficient  of  friction.  Rotation  of  quantities  in 
local  coordinates  into  global  coordinates  is  facilitated  by  the  followiny 
transformations: 
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The  determination  of  the  contactor  code  for  the  present  Iteration 
depends  upon  its  value  for  the  previous  iteration. 

If  i^  was  equal  to  o,  then  node  k was  not  in  contact  during  the  previous 
interation.  In  this  case  we  must  determine  whether  or  not  node  k has  pene- 
trated a target  segment  during  the  present  iteration.  Let  xk  denote  the 
location  of  node  k and  let  , x£,  x£+1  denote  the  locations  of  consecutive 
target  nodes  &-1,  t,  t+1,  respectively,  where  l designates  an  interior  node 
of  some  segment.  We  assume  that  the  entire  list  of  interior  target  nodes 
has  been  searched  and  k is  found  to  be  closest  to  i at  the  end  of  the  present 
iteration.  For  node  i the  interior  of  the  target  is  defined  to  be  that  part 
of  the  plane  consisting  of  the  two  straight  lines  emanating  from  x,  through 

**  ^ i 

I 

?i-l  an<*  -t+l  an<*  ^tending  to  Infinity,  and  all  points  to  the  right  of  1 

these  lines  with  respect  to  the  target  direction.  The  exterior  is  the 
remaining  portion  of  the  plane  (see  Fig.  1-4).  At  the  end  of  an  iteration,  4 

If  x^  is  in  the  Interior  of  the  target  we  say  that  tentative  contact  has 
been  made.  To  determine  if  this  has  occurred  we  employ  the  following 
algorithm  (see  Fig.  1-4  for  notation). 

Let  A = (x^,  » x^,  x^)  and  define  TEST  as  indicated  in  Fig.  1-5. 

If  the  outcome  of  TEST  is  true  (I),  then  x^  is  exterior  to  the  target, 
whereas  if  the  outcome  is  false  (F),  x^  is  in  the  interior  and  tentative 
contact  has  been  made.  In  the  latter  case  further  calculations  are  required 
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to  determine  if  indeed  contact  has  occurred*  and,  if  so,  where.  The 
routine  to  carry  out  these  calculations  requires  the  input  parameters i|nax 
and  which  are  the  maximum  number  of  iterations  allowed  to  determine 

the  approximate  location  of  a contact  point,  and  the  maximum  number  of 
changes  of  target  reference  node  allowed,  respectively.  If  tentative 
contact  has  been  made  we  next  employ  a binary  search  procedure  to  determine 
a good  approximation  to  the  configuration  ft  at  which  initial  contact  was 
actually  made.  With  this  configuration  determined,  we  ascertain  whether 
node  k actually  contacted  the  segment  defined  by  target  nodes  i-1,  i,  t+1 
( 1 a | < 1)  or  did  not  (jaj  > 1).  In  the  former  case  we  set  the  new  value 
of  the  contactor  code  i^  to  1,  whereas  in  the  latter  we  may  try  another 
reference  point  and  repeat  the  calculation.  However,  if  j has  been 
reached,  or  the  new  reference  point  is  a target  segment  boundary  node,  we 
assume  no  contact  has  been  made  and  set  ik  to  o. 

The  following  is  a brief  description  of  the  main  points  of  the  flowchart 
depicted  in  Figs.  1-6  to  1-11. 

If  i^  was  greater  than  zero  for  the  previous  iteration,  then  contactor 
node  k was  in  contact.  In  this  case  it  is  first  checked  if  the  updated 
normal  component  of  traction  xn  is  compressive.  If  this  is  not  the  case, 
then  i^  is  set  to  zero.  If  xn  is  compressive,  then  it  is  determined 
whether  contactor  node  k was  sticking  = 1)  or  whether  it  was  sliding 
(ik  = 2)  during  the  last  iteration.  If  it  was  sticking  then  the  updated 
value  of  | t | is  compared  with  x .t  = f jx_|,  where  f is  the  static 
coefficient  of  friction.  If  |ts|  * Ccrit»  ^en  ^ 1S  se^  eclua1  to  2* 
otherwise  ib  is  set  to  1.  If  node  k was  slidinq  durina  the  last  iteration 
then  a new  value  of  a is  computed  based  upon  the  updated  configuration.  If 

★ 

A contactor  node  can  enter  the  interior  without  passing  through  the  target, 

e.g.  by  sneaking  around  a boundary  node. 
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| a j s 1,  then  | t$ | is  compared  with  Tcrjt  = |-rn  | . The  result  of  this 
comparison  determines  whether  contactor  node  k is  sliding  or  sticking,  as 
above.  If  |ct | •>  1 , the  target  reference  node  is  changed  to  the  appropriate 
adjacent  one.  A maximum  of  icmav  changes  of  reference  are  allowed.  (The 
default  value  of  i’  is  one).  If  icma„  is  exceeded,  the  computation  is 
terminated  and  an  error  message  is  printed.  If  in  changing  the  target 
reference  node  a target  segment  boundary  node  is  encountered,  the  no  contact 
codei^  is  set  to  zero.  If  a value  of  a is  found  such  that  |a|  s 1,  then 
| t$ | is  compared  with  t -t  = f^  |tn|  and  we  proceed  as  described  above. 

When  contactor  nodes  slide  over  target  nodes,  the  computation  of  the  new 
values  of  a is  approximate,  unless  the  target  segment  is  nat.  Thus  some 
caution  is  advised  in  application  to  problems  in  which  substantial  sliding 
is  likely  during  a time  (load)  step. 
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4.  Impact-Release  Conditions 

Imposing  correct  impact  and  release  conditions  are  essential  ingredients 
in  the  accurate  solution  of  dynamic  contact  problems.  The  present  developments 
extend  earlier  work  (see  [1-5])  in  which  simple  node  on  node  normal 
incidence  was  considered.  In  this  section  we  describe  procedures  for  setting 
impact-release  conditions  for  the  stick,  slip  and  frictional  cases.  Since 
each  of  the  cases  differs  somewhat  from  the  others  we  discuss  them  one  at  a 
time.  The  stick  case  is  the  most  straightforward  and  we  shall  describe  it 
first. 

a.  Stick  Contact  Condition 

2 

We  consider  the  case  of  an  open  target  segment  consisting  of  N nodes 
( <e  superscript  refers  to  body  number  2).  The  case  of  a closed  target 
follows  trivially.  We  allow  for  the  possibility  of  an  arbitrary  number 
of  contactor  nodes  impacting  and/or  releasing  the  target  segment  over  the 
time  step.  The  updating  of  nodal  velocities  is  determined  from  the 
following  two  conditions: 

(1)  The  velocity  of  a contactor  node  in  contact  with  the  target 
segment  is  the  linear  interpolate  of  the  target  node  velocities  of  the 
element  boundary  in  question.  For  instance,  consider  the  configuration  of 
Fig.  1-2  in  which  contactor  node  k is  in  contact  with  the  element  boundary 
defined  by  target  nodes  i and  £+1.  In  this  section  we  will  attach  a sub- 
script to  the  nondimensional  location  parameter  to  indicate  that  it  is 
associated  with  contactor  node  k,  i.e.  we  denote  it  c^.  If  u^  denotes  the 
velocity  vector  of  contactor  node  k,  and  u£  and  u£+^  are  the  velocity  vectors 
of  target  nodes  £ and  £+1 , respectively,  then  the  condition  of  linosr 
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interpolation  requires 
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Observe  that  this  definition  remains  meaningful  when  there  are  more  than  one 
contactor  nodes  In  contact  with  a single  element  boundary. 

(2)  The  second  condition  is  a discrete  impulse-momentum  balance. 

Let  A denote  the  tributary  area  surrounding  target  node  i.  Specifically, 

JC 

this  area  is  defined  to  be  one-half  the  length  of  the  two  adjacent  target 

element  boundaries  to  node  i,  in  case  node  i is  an  Interior  node;  and  one- 

half  the  length  of  the  adjacent  target  segment,  if  node  i is  a target  segment 

boundary  node  (see  Fig.  1-12).  Let  denote  the  set  of  contactor  nodes  in 

contact  with  A£  at  the  end  of  the  time  step  (last  iteration),  and  let  C'1 

denote  the  set  of  contactor  nodes  which  were  in  contact  with  A^  at  the 

beginning  of  the  time  step,  but  which  released  during  the  time  step.  Then 

2 

for  each  i € {1,2,...,N  } we  require  that 


u'1 
l ~a 


+ z 

k€c; 


M°  u"1  + M. 
"k  ¥ k 2 


r -i ! 
k€  Ct 


= Mj  u+  + i 
i -i 


« • 


k€C. 


and  for  each  k€C~  we  require 


yC  •+  _ yC  *-1  /it  "I 

Mk  -k  ' Hk  - k ' 2 Ik  ’ 


(1-5) 


U-6) 


where  the  superscripts  + and  -1  indicate  the  updated  values  and  values  from 

4- 

the  previous  time  step,  respectively,  indicates  the  lumped  mass  coefficient 
of  target  node  i and  m£  is  the  lumped  mass  coefficient  of  contactor  node  k. 

Equation  (I-G)  defines  u£  for  all  k€  C~\  In  terms  of  the  data  from 
the  previous  step;  namely  u’^  and  t~^.  Equations  (6.1)  and  (6.2)  lead  to 
the  following  system  of  equations  for  the  target  nodal  velocities: 


A U = B 


(1-7) 


where  A is  tridiagonal  with  nonzero  coefficients 


a.  „ , = r M?(l-aJ,  £t{2,3 N2} 

£»  k£  fc‘i  * k 

kt  ^Vleft 


3£,£  = Mt  + Z 


Ok  + Z 


k^  ^CPleft 


<CI>r1ght 


MjO-ak),  *€  (1,2,3,. ..IT) 


£,  £+1 


k*  bright 


m£  ak,  «•€  (1,2,3,. 


He 


N2-1} 


bn  bi2 


b21  b22 


bN2l  bN22 


in  which 

b ^mJuT1  + z u”1  + At  E x"1,  £6  {1,2,...,  N2} 

£ £ -£  k€  c;  k ~k  2 k€  c;  ~k 

The  subscripts  'left'  and  'right'  on  the  C~‘s  indicate  the  subsets  of 
C to  the  left  and  right  of  target  node  £,  respectively. 

The  second  subscript  on  the  entries  of  U and  B refer  to  the  coordinate 
direction,  e.g.  b^  is  the  Z2  - component  of  b^. 

An  argument  which  employs  dynamic  force  balances  in  place  of  impulse- 
momentum  conditions  yields  relations  for  updated  accelerations  and  contact 
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forces.  The  end  results  are  summarized  as  follows: 


..+ 


• The  updated  accelerations  of  the  target  nodes  are  determined  by 
solving 


A U = B* 


(1-8) 


where 


U = 


J11 


..  + 
U21 


uN21 


12  ; 


••  + i 
U22  ! 


**  o 1 

V2J 


»il 


b21 


B'  = i 


bli2l 


bi2 


b22 


bN22 


i 


b~i  Ki  h k|  c-  Mk  ^k  ’ k|  C-1  Ik 


In  the  above  definition  of  b^,  the  superscripts  refer  to  tne  last  iteration. 

• The  updated  accelerations  of  the  contactor  nodes  which  are  in  contact 
with  the  target  segment  are  given  by  linear  interpolation,  i.e. 


u*.  = + a.  U. 


-k 


V k ~£+l 


• The  updated  accelerations  of  contactor  nodes  which  have  released  from 
the  target  sometime  during  the  time  step  are  given  by 


i - s - ik  / < 


(1-9) 
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• The  updated  contact  forces  associated  with  nodes  which  are  in  contact 
with  the  target  are  given  by 


+ 

3k 


(1-10) 


The  solutions  of  the  systems  of  equations  (1-7)  and  (1-8)  are 
accomplished  most  efficiently  with  an  unsynwietric  tridiagonal  solver.  The 
following  is  a FORTRAN  subroutine  to  carry  out  this  procedure.  The  matrix 
A is  stored  by  rows  as  a one-diinensional  array. 


SUBROUTINE  TRISOL  (A.B.NEQ) 

DIMENSION  A(1),B(NEQ,2) 

C •••  REDUCE  EQUATIONS  TO  UPPER  TRIANGULAR  FORM 
NM  = 3 * NEO-2 
I = 2 

DO  100  N = 1,  NM,  3 
IF  (A(N).  EQ.  0.0)  GO  TO  100 
AA  = A(N+2)/A(N) 

A(N+3)  = A(N+3)  - A(N+1 ) * AA 
B(I,1)  = B(I,1)  - AA  * B ( I - 1 , 1 ) 

B( I ,2)  = B(I,2)  - AA  * B(I-1,2) 

100  1 = 1 + 1 
C ...  BACKSUBSTITUTE 
I = NEQ 

200  IF  (A(NM) . EQ.0.0)  GO  TO  210 
B(I,1)  = B ( I , 1 )/A(NM) 

B ( I , 2 ) = B(I,2)/A(NM) 

210  1 = 1-1 

IF  (I.LE.O)  RETURN 
NK  = NM-3 

B( I ,1 ) = B ( I » 1 ) - A(NM+1 ) * B(I+1 , 1) 

B ( I , 2 ) = B( I ,2)  - A(NM+1 ) * B(I+1,  2) 

GO  TO  200 
END 


b.  Sliding  Contact  Condition 

The  impact  and  release  conditions  for  the  sliding  contact  (frictioi  less) 
case  are  similar  to  the  stick  case,  but  only  involve  normal  components. 

To  describe  the  procedure  employed,  we  will  need  to  introduce  some  new 
terminology.  Local  boundary  coordinates  are  normal  (n)  - tangential  (s) 
coordinates  attached  to  each  target  element  boundary  (see  Fig.  1-3). 
Pseudo-normal  coordinatf"  for  an  interior  target  node  l are  normal  (n)  - 
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tangential  (s)  coordinates  with  respect  to  the  line  joining  target  nodes 
£-1  and  £+1,  whereas  for  target  boundary  nodes  they  are  the  same  as  the 
local  boundary  coordinates  (see  Fig.  1-13).  All  computations  of  Impact  and 
release  data  are  done  with  respect  to  the  configuration  determined  by  the 
last  iteration  within  the  time  step. 

The  velocity  update  is  achieved  as  follows: 

Transform  all  u's  and  t's  appearing  ir.  '1-5)  and  (1-6)  to  pseudo-normal 

2 

coordinates.  Solve  (1-7)  where  B is  replaced  by  the  analogous  N xl  vector 
of  normal  (ri)  components.  U will  then  be  the  updated  n-components  of 
velocity  of  the  target  nodes.  (The  s components  are  unaffected  by  this 
process.)  Rotate  the  n,  s-components  into  local  boundary  coordinates.  Obtain 
the  boundary  normal  (n)  components  of  the  contactor  velocities  by  linear 
interpolation  of  the  updated  n-component  target  velocities.  (The  s-components 
are  unaffected.)  For  nodes  which  have  released,  equation  (1-6)  is  to  be 
applied  with  the  n-components. 

Accelerations  and  contact  forces  are  updated  as  follows: 

Rotate  the  u's  and  j's  appearing  in  the  definition  of  B',  (1-9)  and 
(1-10)  into  pseudo-normal  coordinates.  Solve  (1-8)  where  B'  consists  only 
of  the  n-components;  U will  be  the  updated  n-components  of  the  target  node 
accelerations. 

The  s-components  are  unaffected.  Rotate  the  n , s-components  into  local 
boundary  coordinates  and  linearly  interpolate  to  determine  ii£.  The  s- 
components  are  unaffected.  Updated  n-components  of  contactor  nodes  which 
have  released  are  given  by  (1-9)  with  k replaced  by  ri.  Updated  n-component 
of  contact  force  are  given  by  (1-10)  with  k replaced  by  n.  As  before,  the 
contact  force  for  released  nodes  is  zero, 
c . Frictional  Contact  Condition 

The  impact  ano  release  conditions  for  the  frictional  case  are  identical. 


for  the  n-components,  to  the  frictionless  case.  For  the  s-components  some 
modifications  are  required.  In  the  ensuing  discussion  we  consider  only  the 
s-components.  Again  systems  similar  to  (1-7)  and  (1-8)  are  constructed. 

Here,  let  C”  denote  the  subset  of  nodes  in  contact  which  are  not  sliding  at 
the  last  iteration  of  the  time  step,  and  let  C~^  indicate  the  released  nodes 
(including  sliding  nodes).  Then,  in  the  updating  of  velocities,  s-components 
of  should  be  replaced  by  the  s-components  of  <(xj^  + xj^),  where  k is  a 
shear  correction  factor  which  can  be  adjusted  to  accurately  capture  shear 
wave  phenomena.  In  the  present  work  we  assume  for  simplicity  < = 1.  (These 
remarks  pertain  to  equations  (1-5)  and  (1-6).)  The  s-components  of  velocity 
for  the  nodes  in  C”  are  computed  by  linear  interpolation. 

The  updating  of  the  s-components  of  acceleration  and  contact  force 
proceeds  as  follows: 

Formulate  the  s-component  of  B'  using  u^  and  x~.  Then  U will  be  the 
updated  s-components  of  acceleration  for  the  sticking  nodes.  Released  node 
s-accelerations  are  given  by  (1-9),  and  non-sliding  node  updated  s-contact 
forces  are  given  by  (1-7).  For  nodes  which  are  sliding  there  is  no  update 
of  s-components  of  u and  x. 
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CONTACTOR  NODE 


FIG.  1-1  TYPICAL  CONTACTOR  NODE  AND 
TARGET  ELEMENT  BOUNDARY 
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FIG.  1-2  NONDIMENSIONAL  LOCATION 


PARAMETER 


FIG.  1-5  ALGORITHM  TO  DETERMINE  IF  A 
CONTACTOR  NODE  IS  INTERIOR  OR 
EXTERIOR  TO  THE  TARGET 


FIG.  1-13  PSEUDO-NORMAL  COORDINATES  FOR 
AN  OPEN  TARGET  SEGMENT 
CONSISTING  OF  FOUR  NODES 
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1 1 . New  Developments  with  Particular  Reference  to  the  Hertzian  Algorithm. 

1 . Introduction 

In  this  chapter  we  describe  some  improvements  to  the  Hertzian 
algorithm  documented  in  our  previous  works;  see  [1-5]. 

2.  Automatic  Time  Stepping 

Often  in  impact  problems  the  significant  time  scale  during  contact 
is  much  smaller  than  that  when  the  bodies  are  not  in  contact.  An  example 
is  illustrative.  Consider  the  configuration  of  Fig.  II-l.  A simple 
frame  structure,  perhaps  excited  through  ground  motion,  is  v'brating  in 
its  fundamental  mode.  Let  us  take  the  period  of  the  fundamental  mode 
to  be  order  1.  If  in  the  course  of  the  motion  of  the  frame  structure  it 
impacts  the  rigid  wall  adjacent  to  it,  perhaps  representing  a more  massive 
structure,  the  characteristic  time  scale  while  in  contact  will  be  the 
transit  time  through  the  horizontal  member.  This  could  be  orders  of 
magnitude  less  than  the  period  of  the  fundamental  mode  of  the  frame, 
lo  capture  this  phenomenon  a very  small  time  step  would  have  to  be  taken 
compared  with  the  period  of  the  fundamental  frame  mode.  On  the  other 
hand,  a time  step  this  small  would  be  unnecessary  and  inefficient  while 
the  frame  is  not  in  contact. 

To  effectively  accommodate  situations  such  as  the  one  just  described 
an  automatic  time  step  feature  has  been  programmed  in  FEAP.  Three 
different  time  steps  are  read  in  as  input  data.  The  largest  is  employed 
if  no  contact  is  taking  place.  If  contact  occurs  the  intermediate  time 
step  is  employed  and  the  computation  is  repeated.  The  intermediate  r*:p 
is  used  subsequently  until  contact  is  again  made  at  which  time  the 
computation  is  repeated  with  the  smallest  time  step.  The  smallest  time 
step  is  employed  thereafter  throughout  the  contact  phase.  If  the  bodies 
release  the  largest  time  step  is  again  employed,  and  so  on. 
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This  feature  should  allow  us  to  solve  impact  problems  store 


economically  in  the  future. 
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3 . Higher-Order  Impact-Release  Conditions  ? 

In  this  sectiv*  we  shall  describe  higher-order  impact  and  release 
conditions  W the  Hertzian  algorithm.  The  necessity  of  developing  a 
theory  along  these  lines  was  first  alluded  to  in  [3],  Section  I-3-b,  and 
subsequently  in  [4]  and  [5].  The  theory  is  aimed  at  more  accurately 
capturing  the  post-impact  and  post-release  velocity  states;  the  updated 
accelerations  and  tractions  are  computed  as  was  done  previously  (see  [4] 
or  [5]).  The  reason  for  attempting  to  improve  the  velocity  calculations 
is  that  no  account  is  taken  of  the  acceleration  of  the  nodes  in  question. 

This  can  occasionally  lead  to  poor  results  (see  Section  I-3-b  of  [3]). 

The  theory  presented  herein  accounts  for  acceleration  and  represents  a 
negligible  amount  of  additional  computational  effort. 

We  begin  by  quoting  the  discrete  impact  and  release  conditions  for 
a typical  pair  of  candidate  contact  nodes,  presented  in  [4]  and  [5]: 
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The  notation  is  as  follows:  Superscripts  refer  to  the  body  nuriber  and 
^ = 1,2.  The  subscript  -1  indicates  that  the  quantity  in  question  is 
evaluate  at  the  previous  tine  step,  the  subscript  - refers  to  the  last 
iteration  of  the  present  tine  step,  and  the  subscript  + refers  to  the 
updated  values  accounting  for  impact-release  effects.  Arguments  leading 
to  these  equations  are  presented  in  [1],  [4]  and  [5].  The  refinements 
to  these  quations  to  follow  effecr.  only  the  velocity  equations.  First 
we  consider  the  case  of  impact. 

It  can  be  argued  from  wave-propagation  theory  that  the  updated 
velocity  u+  is  a good  approximation  to  the  velocity  of  the  coalesced 
contact  nodes  at  the  instant  following  impact.  If  impact  occurs  towards 
the  beginning  of  a time  step,  u+  may  not  be  a very  accurate  representa- 
tion for  the  end  of  the  step.  We  seek  to  account  for  this  effect  in  a 
rational  way.  To  do  this  we  make  use  of  the  fact  that  from  the  instant 
after  impact  to  the  end  of  the  time  step,  the  velocity  is  a reasonably 
smooth  function.  With  this  we  define  a new  updated  velocity 

u++  = G+  + u+  At 

where  Jt  = tn+l  - V ^c^^n*  tn+l^  ls  the  instant  of  contact,  and  tp 
and  t^  denote  the  beginning  and  end  of  the  time  step,  respectively. 

The  picture  is  as  illustrated  in  Fig.  I I -2 . It  remains  to  obtain  an 
expression  for  At.  To  do  this,  we  approximate  the  position  of  the 
contact  nodes  in  terms  of  the  data  from  the  previous  time  step,  viz. 

xa  = Xa  + u_“  + (At  - At)  u_“  + i u_“  , 

where  Xa  denotes  the  initial  position  of  the  candidate  contact  mode  of 
body  a.  The  contact  location  is  defined  by 


which  leads  to 


At  = At 


-b  t (b2  - 4ac)1/2 
2a 


where 


“2  "1 
a = u_.|  - u_.j 


b = 


• 2 *1 
U-1  - U-1 


c = X 


X1  ♦ uj 


- u 
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The  physically  relevant  solution  is  determined  by  the  condition 


0 <_  At  £ At 

If  no  solution  satisfies  this  condition  At  is  set  to  At/2. 

To  obtain  improved  release  velocities  we  also  attempt  to  estimate  a 
more  accurate  time  of  release  within  the  step.  Since  the  nodes  were  in 
contact  at  the  end  of  the  previous  step,  t_-j  >0.  Let  t denote  the  last 
value  of  contact  force  before  release  occurred.  In  keeping  with  our 
previous  conventions  (see  [4]  or  [5]),  t will  be  negative  (indicating 
tension)  or  less  than  or  equal  to  2%  of  t_-|.  In  the  latter  case  we 
maintain  the  use  of  u“  as  the  post-release  velocity.  (We  note  that  the 
second  term  on  the  right-hand  side  of  the  expression  for  u+  represents 
the  impulse  over  the  time  step,  assuming  linear  interpolation  between 
t_i  and  zero.)  In  the  case  in  which  T < 0,  we  compute  (see  Fig.  1 1 -3 ) 

tv  - z nz  . \ 

“ l/\t  — t -t/ 


The  new  updated  velocity  is  then  defined  to  be 


III.  Sample  Problem 

1 . Identification  of  Urethane-Polystyrene  Composite  Foam  for  the 
Nonlinear  Continuum  Element. 

The  theory  for  a nonlinear  elastic  continuum  element  was  presented 
in  [3],  Sections  I 1-2  and  I 1-3.  This  element  has  been  programmed  in 
FEAP  and  the  following  options  are  available: 


2-0 


plane  strain 
plane  sti*ess 
axisyimetric 


3-D 


Two  different  quadrature  rules  can  be  employed  in  each  case:  2x2 

Gaussian  quadrature  for  all  terms,  or  2 x 2 Gaussian  quadrature  for 
u-terms  and  1-point  quadrature  for  X-terms  (see  [3]  for  notation  and 
further  details).  The  latter  option  is  appropriate  for  incompressible 
and  nearly-incompressible  applications.  For  use  in  subsequent  check 
problems,  we  have  selected  the  parameters  x and  y so  that  a close  fit 
is  obtained  to  the  loading  cycle  for  a urethane-polystyrene  composite 
foam  described  in  [6].  Values  of  E = 50  psi  and  v = .1,  .25,  .3  were 
selected  which,  for  the  configuration  illustrated  in  Fig.  III-l,  leads 
to  the  response  illustrated  in  Fig.  1 1 1-2 . As  can  be  seen,  the  best 
results  are  obtained  with  v = .25.  In  subsequent  calculations, 
unless  otherwise  indicated,  E = 50  psi  and  y = .25  were  employed. 
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2 . Equivalence  of  Present  Incompressible  Formulation  with  a Mean-Pressure- 
Variable  Element. 

In  [3],  on  p.  32,  we  conjectured  that  employing  2x2  Gaussian  quad- 
rature on  the  v- terms,  and  1 -point  quadrature  on  the  A -terms,  for  the 
standard  four-node  isoparametric  quadrilateral,  might  yield  result"  identical 
to  the  constant  mean  pressure-bilinear  displacement  element  employed  in 
the  past  by  Hughes  and  Allik  [7].  We  have  attempted  to  corroborate  this 
conjecture  by  performing  an  analysis  using  both  elements.  Consider  the 
configuration  illustrated  in  Fig.  III-3.  The  beam  is  modelled  with  plane 
strain  rectangular  using  several  different  quadrature  rules  and  the 
constant  mean  pressure-bilinear  displacement  element  of  [7], 

The  beam  is  fixed  at  the  left  end  and  a uniform  shear  is  applied  to 

the  right  end.  The  results  confirm  the  equivalence  of  the  'underintegrated' 

element  with  the  constant  mean  pressure-bilinear  displacement  element. 

We  thank  H.  Allik  and  P.  Caccistore  of  the  Electric  Boat  Division  of 
General  Dynamics,  Groton,  Connecticut,  for  providing  us  with  the  results 
for  the  constant  mean  pressure-bilinear  displacement  element.  (Recently 
an  analytical  study  has  been  performed  which  establishes  the  equivalence 
of  the  two  elements;  see  Hughes  [8].) 
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3.  Quasi -Static  Analysis  of  a Skull -Pad  Contact  Configuration 

This  analysis  consists  of  a cylindrical  'skull'  model  being  contacted 
with  a soft  pad.  The  material  properties  of  the  three-layer  skull  are 
those  quoted  in  [3],  Table  1-1*,  and  originally  obtained  from  [9].  The 
pad  material  is  that  described  in  Section  III-l.  Linear  elements  are 
used  to  model  the  skull  and  nonlinear  elements  are  used  to  model  the 
pad.  All  elements  employ  the  plane  strain  option  and  2x2  Gaussian 
quadrature.  The  skull  is  fixed  at  the  uppermost  node  and  the  pad  is 
driven  into  the  skull  and  withdrawn  by  way  of  prescribing  a uniform 
displacement  condition  along  the  bottom  of  the  pad.  The  initial  gap 
between  skull  and  pad  is  0.1  inches.  The  maximum  vertical  displacement 
of  the  bottom  of  the  pad  is  0.5 inches  and  the  displacement  is  applied 
in  steps  of  0.1  inches.  Unloading  is  performed  similarly.  The  contact 
condition  is  assumed  to  be  perfect  friction  along  the  contact  surface. 

Thus  there  is  no  tangential  slipping  while  in  contact.  Release  occurs 
when  tension  is  sensed  normal  to  the  target  segment  (in  this  case  the 
pad).  The  analysis  employs  the  'between  node'  contact  element  described 
in  [2],  Section  I 1-3,  and  the  kinematically  nonlinear  search  algorithm 
described  in  [3],  Section  III.  The  target  segment  consists  of  the  seven 
element  boundaries  along  the  top  surface  of  the  pad.  There  is  a total 
of  seven  contactor  nodes  — the  innermost  seven  nodes  along  the  boctom 
outer  surface  of  the  skull.  Initial  and  deformed  configurations  are 
depicted  in  Fig.  II 1-4.  The  unloading  steps  are  identical  to  the  loading 
steps  (i.e.  step  numbers  6-4,  7=3,  8=2,  9=1,  10=0)  and  thus  not  shown. 

Note  that  the  contactor  nodes  contact  the  target  segment  between  nodes. 

* there  is  a typographical  error  in  the  table.  The  density  of  the  brain 
material  should  be  .937  x 10~4  Ib-sec^/in  . This  was  the  value 
actually  used  in  the  analyses. 
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The  strain  in  the  pad  elements  near  the  axis  of  symmetry  reach  a maximum 
strain  of  approximately  502  and  the  pad  elements  near  the  periphery  of 
the  contact  zone  experience  a maximum  rotation  of  approximately  30°. 

Normal  and  tangential  contact  stresses,  at  the  point  of  maximum  contact 
area  development,  are  plotted  in  Fig.  I I 1-5.  The  pad  has  the  effect  of 
more  uniformly  distributing  the  contact  force  than  a rigid  surface;  cf. 

[3],  Section  I-l-a.  Total  vertical  contact  force  and  contact  area  width 
are  plotted,  versus  applied  displacement  step  number,  in  Fig.  I I I -6. 

In  our  initial  attempts  to  solve  this  problem  we  observed  a lack 
of  convergence.  This  was  due  to  the  following  situation:  Nodes  frequently 
released  and  then  recontacted  during  iterating  within  a step.  The  point 
of  contact  was  set  to  the  last  penetration  point,  rather  than  the  initial 
contact  point.  This  was  in  violation  of  the  stipulated  no-slip  condition. 
In  addition,  the  contactor  node  along  the  symmetry  axis  penetrated  the 
contact  surface  without  contact  being  sensed.  This  was  due  to  a small 
negative  horizontal  displacement,  caused  by  round-off,  which  made  the 
search  algorithm  think  the  contactor  node  was  outside  the  target  segment. 
Upon  correcting  these  fallibilities  convergence  was  achieved. 
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u = O.it 


u = O.l  t 


FIG.  II 
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AREA  WIDTH  VERSUS  APPLIED 
DISPLACEMENT 


IV.  A Simple  and  Efficient  Finite  Element  for  Plate  Bending 


1 . Introduction 

An  enormous  amount  of  effort  has  been  devoted  to  the  development  of 
finite  elements  for  the  bending  of  plates.  The  literature  is  extensive 
and  we  will  not  make  an  attempt  to  review  it  here.  (The  interested 
reader  may  consult  any  of  the  standard  texts  for  references,  e.g. 

[10-12].)  Host  of  this  effort  has  been  oriented  towards  linear 

problems;  in  particular,  to  the  classical  Poisson  - Kirchhoff  theory  of 
bending.  The  -continuity  requirement  imposed  by  this  theory  on 
'displacement'  finite  element  models  precludes  the  development  of  simple 
and  natural  elements  (see  [13]).  Because  of  this,  incompatible  elements 

(e.g.  [14-15])  are  oft.-n  resorted  to,  since  they  involve  simpler 
programming  than  the  rather  complicated  C1 -continuous  elements  (e.g. 
[15-17]  ) and  are  competitive  from  an  accuracy  standpoint. 

Accurate  higher-order  (^-elements  have  also  been  developed  (e.g. 
[18-20]),  but  they  too  are  quite  complicated  and  involve 
nodal  derivative  degrees -of- freedom  of  order  greater  than  one,  which 
complicates  the  specification  of  boundary  conditions. 

The  assumed  stress  hybrid  bending  elements  of  Pian  and  his  associates 
(e.g.  [21])  have  proven  to  be  accurate,  but  they  have  some  drawbacks  and 
thus  are  not  widely  used. 

Another  approach  to  the  development  of  bending  elements  for  thin 
plates  involves  the  so-called  ‘discrete  Kirchhoff  hypothesis' 

(e.g.  [22,23]).  In  this  approach  the  classical  equations  are 

abandoned  in  favor  of  a bending  theory  which  includes  shear  deformations . 

The  result  is  that  only  (^-continuity  is  required  of  the  shape  functions. 


To  capture  the  behavior  of  thin  plate  theory,  the  constraint  of  zero  shear 
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strains  is  imposed  at  a discrete  number  of  points.  The  method  is 
effective,  but  implementations  tends  to  be  somewhat  complicated.  Recent 
inprovements  and  variants  on  this  theme  have  been  proposed  by  Fried 
[24-26]. 

An  accurate  quadrilateral  element  for  thick  and  thin  plates  has 
been  developed  by  enkiewicz,  Taylor  and  Too  [27].  This  element  possesses 
eight  nodes--four  corner  and  four  midside--';ith  the  basic  three  degrees- 
of- freedom  per  node.  The  transverse  displacement  and  rotation  shape 
functions  are  selected  from  the  'serendipity'  family  (see  [13])-  Two-by- 
two  Gaussian  quadrature  is  an  essential  requirement  for  the  good  performance 
of  the  element. 

In  summarizing  these  developments  one  can  confidently  assert  that 
for  linear  problems  of  plate  bending  many  adequate  elements  exist.  The 
choice  is  more  a matter  of  taste  as  no  single  element  is  clearly  superior 
to  the  rest  in  all  cases. 

Many  users  of  finite  element  computer  programs  find  a 'basic' 
four-node  quadrilateral  element  particularly  appealing  due  to  its 
simplicity.  It  is  our  feeling  that  this  appeal  will  become  even  greater 
<hen  nonlinear  applications  are  undertaken.  In  the  nonlinear  regime— 
and  especially  in  nonlinear  dynamics— computational  cost  is  the  prime 
concern.  Due  to  frequent  reformulations  of  tangent  stiffnesses,  complicated 
element  routines  can  lead  to  exorbitant  computational  expenditures  and 
may  actually  preclude  nonlinear  analysis.  A simpler  element  of  competi- 
tive accuracy  becomps  quite  desirable  under  such  c i crums tances . Other 
factors  in  nonlinear  analysis  buttress  tnis  assertion.  For  example,  the 
accuracy  level  attainable  in  nonlinear  problems  is  often  severely  limited 
due  to  the  uncertainty  of  nonlinear  material  characterizations.  Thus  it 
makes  little  sense  to  engender  significant  computational  cost  for  complicated 
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bending  elements  which  are  only  marginally  more  accurate  than  simpler 
elements,  since  the  confidence  level  of  the  overall  analysis  may  be 
affected  only  negligibly.  Unfortunately,  heretofore,  no  really  simple 
alternative  has  existed. 

In  this  chapter  we  attempt  to  remedy  this  situation.  We  develop 
what  we  believeis  the  simplest  effective  plate  bending  element  yet  pro- 


pose^. The  element  is  a four-node  quadrilateral  with  the  basic  three 
degrees -of- freedom  per  node.  The  element  shape  functions  are  bilinear 
for  transverse  displacement  and  rotations.  The  shear  'locking' 
associated  with  such  low-order  functions  in  application  to  thin  plates 
is  eleviated  by  splitting  the  shear  and  bending  energies  and  using  one- 
point  quadrature  on  the  shear  term.  The  simplicity  of  the  element  lends 
itself  to  concise  and  efficient  computer  implementation. 

To  develop  the  theory  in  its  simplest  setting,  we  consider  in 
Section  IV-2  a beam  element  involving  linear  displacement  and  rotation 
shape  functions.  We  show  how  exact  integration  (two-point  Gaussian 
quadrature)  of  the  element  stiffness  matrix  leads  to  an  overly  stiff 

: element  and  we  present  an  heuristic  argument  why  this  is  the  case.  We 

then  show  how  employing  one-point  quadrature  on  the  shear  term  lessens 

the  stiffness.  The  concept  is  identical  for  the  plate  bending  element 

i 

' which  is  developed  in  Section  I V- 3 . The  effectiveness  of  the  element 

i 

i 

\ in  thin  plate  bending  is  demonstrated  in  Section  IV-4.  A simple  computing 

strategy  for  dealing  with  the  numerically  sensitive  case  of  extremely 
thin  piates  is  presented  in  Section  I V - 5 . In  Section  IV-6  we  consider 

l applications  to  thick  plates.  It  is  shown  that  the  element  is  still 

f effective  for  moderately  thick  plates.  However,  for  very  thick  plates, 

in  which  the  thickness  of  individual  elements  exceed  their  characteristic 

I 

? lengths,  a slight  modification  of  t.he  shear  quadrature  need  be  employed. 

L— .. 
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2.  Example:  Linear  Beam  Elemant 


In  this  section  we  will  describe  the  formulation  of  a beam  element 
stiffness  for  which  displacement  and  rotation  are  assumed  to  be  indepen- 
dent linear  functions.  Exact  integration  of  the  element  stiffness  can  be 
facilitated  by  two-point  Gaussian  quadrature,  whereas  one-point  integration 
exactly  integrates  the  bending  contribution,  but  'underintegrates'  the 
shear  contribution.  For  the  case  of  thin  beams  we  view  the  shear  term 
as  a constraint  which  attempts  to  enforce  the  condition  of  negligible 
shear  strains.  We  shall  show  that  one-point  quadrature  has  a decisively 
positive  effect  on  the  accuracy  of  the  element;  two-point  quadrature 
leading  to  worthless  numerical  results.  The  upshot  of  all  this  is  that 
by  appropriately  underintegratinq  troublesome  terms,  good  bending  behavior 
can  be  attained  by  the  simplest  shape  functions. 

The  equations  of  a rectangular  cross-section  beam,  including  shear 
deformation  effects,  emanate  from  the  following  expression  for  strain 
energy: 


U(w,e) 


(IV-1) 


where  w is  the  transverse  displacement  of  the  center-line,  0 is  the 
rotation  of  the  cross-section,  E is  Young's  modulus,  G is  the  shear 
modulus,  < is  the  shear  correction  factor  (throughout  we  employ  *c  = 5/6), 
t is  the  depth,  L is  the  length  and  x is  the  axial  coordinate.  The  first 
term  on  the  right-hand  side  of  (IV-1)  is  the  bending  energy  and  the 
second  is  the  shear  energy.  With  independent  expansions  for  w and  e, 
(IV-1)  can  be  employed  to  derive  beam  element  stiffness  matrices.  The 
case  we  are  interested  inis  when  both  w and  e are  assumed  to  behave 
linearly  over  an  element.  This  leads  to  a four-degree-of- freedom  element 
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in  which  displacement  and  rotation  are  the  nodal  degrees-of-freedom.  By 
virtue  of  the  fact  de/dx  is  constant  within  this  element,  the  bending 
energy  may  be  exactly  evaluated  by  one-point  Gaussian  quadrature.  On 
the  other  hand,  two-point  Gaussian  quadrature  is  required  to  exactly 
integrate  the  shear  energy  term  due  to  the  explicit  presence  of  e,  which 
is  linear  within  the  element.  Employing  one-point  quadrature  on  the 
shear  energy  term  ’underintegrates'  the  element  and  it  is  our  prime 
concern  here  to  ascertain  the  effect  of  this  procedure.  (See  also 
Gallagher  [11],  pp.  364-367.) 

A series  of  test  computations  were  performed  to  determine  the 
behavio/  of  the  element.  A cantilever  beam  subjected  to  an  end  load 
(see  Fig.  IV-1)  was  analyzed  for  various  discretizations.  The  first 
example  is  for  a relatively  deep  beam.  The  data  are: 

E = 1000 
G = 375 

t = 1 

L = 4 

Tip  displacement  results  for  several  discretizations  are  presented  in 
Table  IV-1.  As  is  evident,  the  one-point  quadrature  results  are  vastly 
superior  to  the  two-point  results.  A more  severe  test  for  linear 
elements  is  bending  governed  by  Bernoulli  - Euler  theory.  In  this  case 
shear  strains  are  to  be  equal  to  zero.  Such  a situation  can  be  brought 
about  in  the  present  theory  if  depth  is  taken  very  small  compared  with 
element  length.  Alternatively,  a very  large  fictitious  value  of  G can 
be  specified.  In  the  second  example  we  attempt  to  ascertain  the  behavior 
of  the  linear  element  when  the  assumptions  of  the  Bernoulli  - Euler 
theory  apply.  The  data  of  the  previous  example  are  employed  with  the 
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Table  IV-1.  Normalized  tip  displacement 
for  a deep  cantilever  beam. 


Tip  displacement— 
one-point  quadrature 


Tip  displacement- 
two- point  quadrature 


.762 

.416  x 10' 

.940 

.445 

.985 

.762 

.996 

.927 

.999 

.981 

Table  IV-2.  Normalized  tip  displacement 
for  a thin  cantilever  beam. 


Tip  displacement— 
one-point  quadrature 


Tip  displacement— 
two-point  quadrature 


.2  )0  x 10 
.800  x 10' 
.320  x 10' 
.128  x 10' 


.312  x 10' 
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exception  of  G which  is  set  here  to  375  x 10  . Results  are  listed  in 
Table  IV-2.  The  one-point  quadrature  results  are  quite  accurate  whereas 
the  two-point  results  are  in  error  by  approximately  three  orders  of 
magnitude.  Early  attempts  at  developing  bending  elements  with  simple 
shape  functions  were  abandoned  because  of  results  like  those  for  the 
two-point  quadrature  presented  here. 

We  shall  now  proceed  to  give  a heuristic  argument  why  two-point 
quadrature  causes  such  an  overly  stiff  element.  Consider  a cantilever  beam 
discretized  into  N elements.  In  the  assembled  stiffness  matrix  there  are  2N 
Jegrees-of- freedom  --  two  degrees-of- freedom  per  element.  The  shear  contribution 
to  the  stiffness  represents  a constraint  on  the  shear  strains.  If  one- 
point  quadrature  is  employed,  one  constraint  is  imposed  upon  the  element, 
whereas  if  two-point  quadrature  is  employed,  two  constraints  are  imposed 
upon  the  element.  In  the  latter  case  the  number  of  constraints  per 
element  equals  the  number  of  degrees -of- freedom  per  element,  and  the 
result  is  that  the  mesh  'locks'. 

This  can  be  seen  more  precisely  by  looking  at  the  stiffness 
contributions  of  the  bending  and  shear  terms.  We  assume  the  nodal 
degrees  of  freedom  are  ordered  as  follows:  w-| , 6p  *2>  °2»  anc*  ^ 

the  element  length.  The  stiffnesses  are: 
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(IV-2b) 


(IV-2c) 


where  k^  is  the  bending  stiffness,  and  and  k^  are  the  one-point 

and  two-point  quadrature  shear  stiffnesses,  respectively.  It  is  easily 
verified  that  the  rank  of  k^  is  one  and  the  rank  of  kj^  is  two.  In 
the  latter  case,  k^  is  completely  dominated  by  the  shear  stiffness,  as 
the  following  simple  example  illustrates. 

Consider  the  case  of  a one-element  cantilever  beam,  subjected  to  an 
end  load  P and  moment  M. 

( i ) One-point  quadrature. 

Combining  (IV-2a)  and  (IV-2b),  eliminating  appropriate  rows  and 
columns,  and  solving  for  the  tip  displacement  and  rotation  yields 


w = (h2/4c-  + ff 1 ) P + hM/2a  , 


0 = (hP/2  + M)/a  , 
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( I V- 3a ) 
( I V- 3b ) 
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where 

a = Et3/12h, 
g = *Gt/h . 

In  the  thin  beam  limit  (i.e.,  when  g»a),  (IV-3a)  becomes 

w = h(hP/2  + M)/2o, 


(IV- 3c) 
(IV-3d) 


(IV-4) 


and  ( I V - 3b ) remains  unchanged.  Thus  we  are  left  solely  with  the  deforma- 
tion due  to  bending  as  is  right. 

(ii)  Two-point  quadrature. 

Carrying  out  the  same  steps  as  in  case  i,  with  (IV-2c)  in 
place  of  (IV-2b),  yields 

.2, 


where 


w = P + hM/2y, 


G = (hP/2  + M)/y, 


y = a + h‘g/12. 


(IV-5a) 

(IV-5b) 

(IV-5c) 


In  the  thin  beam  limit  (IV-5a)  and  (IV-5b)  become 


w = (4P  + 6M/h)/B  v 


o = 6 ( hP  + 2M)/h^B , 


(IV-oa) 
( I V-6b) 


l 


respectively.  In  this  case  only  deformations  due  to  shear  are  in  evidence 
and  (IV-6a)  and  ( I V-6b ) are  Q(t  3)  in  error. 

In  passing  we  note  that  there  are  some  circumstances  in  which  the 
present  element  may  have  some  practical  value.  For  example,  an  axisymme- 
tric  shell  version  miqht  be  useful  for  shell  covered  solids  in  which 
bilinear  elements  are  used  to  model  the  solid.  The  fact  that  only  one 


The  strain  energy  for  an  isotropic,  linearly  elastic  plate,  including 
shear  deformation,  is 

U(w , > 89 ) - 


where  and  ^ are  cartesian  coordinates,  w is  the  transverse  displacenent, 
0 1 and  are  the  rotations  about  the  x-|  and  x^  axes,  respectively,  E is 
Young's  modulus,  v is  Poisson's  ratio,  < is  the  shear  correction  factor, 
t is  the  plate  thickness  and  A is  its  area.  The  first  integral  in  (IV-7) 
represents  the  bending  energy  and  the  second  represents  the  shear  energy. 

'.ie  consider  a four-node  quadrilateral  element  and  assume  the  displacement 
and  rotations  are  expanded  in  independent  bilinear  shape  functions.  The 
isoparametric  concept  is  employed  (see  Zienkiewicz  [10]).  This  results 
in  three  degrees -of -freedom  --  one  displacement  and  two  rotations  --  at  each 
of  the  comers. 

For  very  thick  plates  two-by-  jho  Gaussian  quadrature  leads  to 
acceptable  results,  however,  for  thin  plates  it  causes  'locking'  as 
indicated  for  the  beam  in  the  previous  section.  In  this  case  we  use 
two-by-two  Gaussian  quadrature  on  the  bending  energy  term  and  one-point 
iaussian  quadrature  on  the  shear  energy  term.  This  results  in  one  con- 
straint per  element.  In  large  meshes  there  are  approximately  three 
equations  per  element,  tnus  there  is  no  danger  of  the  mesn  'i^cring'. 


As  is  apparant,  cbe  proposed  element  is  extremely  simple,  and  easily  and 
concisely  coded.  We  are  certain  that  the  element  routines  are  faster 
than  any  other  plate  bending  element  yet  proposed.  In  the  next  section 
we  will  show  that  the  element  is  also  surprisingly  accurate. 
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4.  Numerical  Examples:  Thin  Plates 


In  this  section  we  present  several  numerical  examples  which  have 


become  more  or  less  standard  ones  for  evaluating  plate  elements.  All 
computations  were  performed  on  a CDC  6400  computer  in  single  precision. 
(A  single  precision  word  consists  of  60  bits  on  the  CDC  6400.) 


Square  Plate 


The  data  for  this  example  consists  of  the  following  (see  Fig.  IV-2) 


E = 10.92  x 10' 


Both  simply  supported  and  clamped  boundary  conditions  were  considered  as 


well  as  concentrated  and  uniformly  distributed  loadings.  Results  are 


presented  in  Tables  IV-3  and  IV-4  for  < values  of  1000  and  5/6.  The 


former  value  is  set  to  artificially  maintain  the  Poisson- Kirchhoff  con- 


straint. Due  to  the  fact  that  the  plate  is  rather  thin  (L/t  = 100), 


there  is  little  difference  in  the  results  for  the  two  values  of  < , In 


fact,  the  bending  moments  are  identical.  In  practical  situations  there 


seems  no  point  in  exceeding  the  ‘natural’  value  of  < - 5/6. 


Moment  and  shear  resultants,  and  displacements  are  plotted  >n  Fig. 


I V- 3 along  the  line  = 0,  for  the  clamped,  concentrated  load  case  in 

./hich  r = 5/6.  Along  x-j  = 0 arid  x^  = 0,  the  x-j  and  X£  components  of  the 


moment  are  equal,  as  are  the  x ) and  components cf  the  shear. 


The  simply  supported  concentrated  load  case  has,  it  seems,  taken 


on  the  role  of  the  preeminent  comparison  problem  for  bending  elements. 


In  Fig  IV-4  the  present  element,  with  < = 1000,  is  compared  with  data 


taken  from  Gallagher  [!:].  The  good  convergence  of  the  element  is  evident. 


Clamped  Circular  Plate 


The  data  for  this  example  is  the  sane  as  for  the  previous  problem 
except  (see  Fig.  IV-5) : 


R = 5 
t = .1 

Results  are  presented  in  Table  IV-5  for  concentrated  and  uniform  loadings, 
and  * values  of  1000  and  5/6. 

Again,  due  to  the  thinness  of  the  plate,  there  is  little  difference 
in  the  displacement  results  for  the  two  values  of  <,  and  the  moment  results 


l»li] 


Table  IV-3.  Normalized  center  displacement 
and  bending  moment  for  a simply 
supported  square  plate. 


Number  of 
elements 


Number  of 
element 


Displacement- 
Concentrated  load 

Displacement- 
Uniform  load 

Moment- 
Uniform  load 

.9922 

.9770 

.351 

.9948 

.9947 

.963 

.9982 

.9987 

.991 

(a)  K = 1000 

Displacement- 
Concentrated  load 

Displacement- 
Uniform  load 

Moment- 
Uniform  load 

.9957 

.9782 

851 

.9991 

.9960 

.063 

1.0034 

l 

.9997 

.991 

= 5/6 


m 


Table  IV-4.  Normalized  center  displacement 
and  bending  moment  for  a clamped 
square  plate. 


Nunber  of 
elements 


Displacement— 
Concertrated  load 

Displacement- 
Uniform  load 

Moment- 
Uniform  load 

.8652 

.9535 

.822 

.9650 

.9850 

.955 

.9920 

.9937 

.986 

Number  of 
elements 


Displacement- 
Concentrated  load 

Displacement- 
Uniform  load 

Moment- 
Uniform  load 

.8720 

.9575 

.822 

.9748 

.9890 

.955 

1.0034 

.9976 

.986 

Table  IV-5.  Normalized  center  displacement 
and  bending  moment  for  a clamped 
circular  plate. 


Number  of 
dements 

Displacement- 
Concentrated  load 

Displacement- 
Uniform  load 

Moment- 
Uniform  load 

3 

.9197 

.8587 

.827 

12 

.9579 

.9535 

.957 

.9883 

.9888 

.990 

Nuiriber  of 
elements 

Displacement- 
Concentrated  load 

Displacement- 
Uniform  load 

Moment — 
Uniform  load 

3 

.9267 

.8621 

.327 

12 

.9674 

.9570 

.957 

48 

1 .0005 



.9925 

-990 

- - 

. = 5/6 
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S.  Numerical  Sensitivity  due  to  Extreme  Thinness 


The  results  of  the  previous  section  indicate  that  despite  the 
simplicity  of  the  present  element  it  is  quite  accurate.  However,  one 
precaution  must  be  taken  when  employing  elements  derived  from  the  reduced 
integration  concept.  This  admonition  stems  from  the  observation  that 

o 

the  shear  stiffness  is  0 ( (h/t)  ) times  the  bending  stiffness.  (In  the 
case  of  a quadrilateral  bending  element  h may  be  taken  to  be  the  length 
of  the  longest  side.)  For  fixed  h,  as  t -*•  0,  it  is  only  a matter  of  time 
before  the  effect  of  the  bending  stiffness  vanishes  completely  due  to  the 
finite  computer  word  length.  Results  of  this  kind  can  be  seen  for  the 
beam  element  in  Fig.  IV-6  and  for  the  plate  element  in  Fig.  IV-7.  The 
plateaus  represent  the  appropriate  solutions  for  the  meshes  in  question 

in  the  'thin'  limit.  Deterioration  of  the  numerical  solution  begins  to 

4 5 

occur  at  L/t  = 10  in  the  case  of  the  beam,  and  L/t  = 10  for  the  plate; 

this  corresponds  to  element  aspect  ratios  (i.e.,  h/t)  of  1 0^/16  and 

5 

10  /8,  respectively.  It  is  unlikely  that  aspect  ratios  larger  than 
these  values  will  be  met  in  practice.  However,  on  computers  with  shorter 
single  precision  word  length, deterioration  will  commence  at  smaller 
aspect  ratios.  Here  it  is  important  to  employ  a strategy  which  circum- 
vents these  difficulties.  This  can  be  done  as  follows:  Determine  the 

maximum  element  aspect  ratio  for  which  good  results  are  obtained  by 
plotting  graphs  similar  to  Figs.  IV-6  and  IV-7.  Before  combining  the 
shear  and  bending  contributions  of  the  element  stiffness  test  the  aspect 


« r 1 ^ 

Id  I C 3 3 
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fn  iv.ii  yv/vfu  i cau  i v d aft:  v/utu 


combine  in  the  usual  way.  Otherwise,  multiply  the  shear  stiffness  by 
2 

(t/h)  times  the  square  of  the  maximum  element  aspect  ratio  allowed, 
then  confcine.  This  will  reduce  the  disparity  between  the  bending  and 
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shear  term  to  an  acceptable  level,  yet  maintain  thin  element  behavior. 
Numerical  results  illustrating  these  ideas  are  depicted  in  Figs.  IV-6 
and  IV-7.  The  maximum  allowable  aspect  ratios  were  determined  to  be 
10^/16  for  the  beam  element  and  10^/8  for  the  plate  element,  from  Figs. 
IV-6  and  IV-7,  respectively.  Employing  these  values  in  the  procedure 
described  above  enables  the  plateaus  to  be  extended  to  the  higher  aspect 
ratios,  as  illustrated  in  Figs.  IV-6  and  IV-7.  This  process  enables  the 
reduced  integration  procedure  to  be  applied  successfully  in  cases 
involving  arbitrarily  large  aspect  ratios. 

Of  course,  another  way  to  avoid  difficulties  is  to  work  in  double 
precision  on  short  word  computers. 


"vju'aifi 


r 
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6.  Application  to  Thick  Plates 

One-point  Gaussian  quadrature  on  the  shear  term  has  a decidedly 
beneficial  effect  in  application  to  thin  plates.  However,  the  opposite 
is  true  for  very  thick  plates.  The  difficulties  do  not  always  manifest 
themselves.  For  example,  results  for  the  uniformly  loaded  clamped 
circular  plate  are  acceptable  (see  Fig.  IV-8;  q refers  to  the  magnitude 
of  the  load).  For  t/R  = .4  the  results  are  of  about  the  same  level  of 
accuracy  as  those  in  [27],  where  a higher-order  element  is  employed. 

On  the  other  hand,  results  for  the  same  plate  subjected  to  a concentrated 
load  tend  to  oscillate  about  the  exact  solution  (see  Fig.  IV-9).  This 
problem  is  a trying  one  as  the  exact  Reissner's  theory  solution  consists 
of  an  infinite  displacement  under  the  load,  viz. 


bending  shear 


where  P is  the  concentrated  force  and  D is  the  bending  rigidity. 

As  the  plate  thickness  is  reduced  the  oscillations  are  lessened  (see 
Figs.  IV-10  to  I V- 1 5 ) . From  these  results  we  conclude  that  when  the 
t/h  ratio  exceeds  unity,  the  one-point  Gaussian  quadrature  of  the  shear 
term  should  be  abandoned  in  favor  of  the  following  scheme:  Two-by-two 

Gaussian  quadrature  should  be  used  on  the  (aw/ax^)~  and  (aw/ax^) 
contributions  to  the  shear  enerqy.  The  remaining  terms  in  the  shear 
energy  should  be  integrated  as  usual  by  one-point  Gaussian  quadrature. 
We  refer  to  this  element  as  the  ’ r-od i f i ed ' one-point  shear  element. 


L. 


A spectral  analysis  of  the  element  stiffness,  when  one-point 
Gaussian  quadrature  is  employed  on  the  shear  term,  reveals  that  there 
are  five  zero  eigenvalues  --  two  more  than  the  usual  three  rigid  body 
modes.  Thus  the  element  by  itself  is  rank  deficient,  but  this  only 
manifests  itself  in  problems  for  very  thick  plates  and  here  only  in 
certain  cases.  The  two  additional  zero-jnergy  modes  are  illustrated 
in  Fig.  IV-16.  The  first  mode  consists  of  e-j  = e2  = 0 and  an  'hourglass' 
pattern  for  w (see  [38]).  Modifying  the  one-point  shear  integration  as 
indicated  above  removes  the  hourglass  mode  and  leads  to  good  results  for 
very  thick  plates  as  evidenced  by  Figs.  IV-9  and  IV-10.  The  second  mode 
consists  of  w = 0 and  an  in~plane  twisting  of  the  plate.  In  a mesh  in 
which  the  rigid  body  modes  are  removed,  this  pattern  cannot  persist  and 
thus  causes  no  further  rank  deficiency. 
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7.  Conclusions 

We  have  presented  an  element  for  the  bending  of  thin  and  moderately 
thick  plates  which  involves  minimal  programming,  is  highly  efficient  and 
competitively  accurate.  Due  to  these  attributes  the  element  offers  an 
attractive  basis  for  nonlinear  developments.  Numerical  sensitivity  in 
applications  involving  extremely  thin  elements  has  been  shown  to  be 
avoidable  by  employing  a simple  computational  strategy.  Very  thick 
plates  may  be  successfully  analyzed  by  a slight  modification  to  the 
element. 


FIG.  IV-2  SQUARE  PLATE.  DUE  TO 

SYMMETRY  ONLY  ONE  QUADRANT  IS 
DISCRETIZED 


Blti  i 


1.0  2.0  3.0  4.0  5.0 


AND  DISPLACEMENT  (w)  VERSUS  EDGE 
COORDINATE  (wj ) TOR  A CLAMPED 
SQUARE  PLATE  SUBJECTED  TO  A 
CONCENTRATED  LOAD 
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LOAD.  COMPARISON  Or  CHUR 
DISPLACEMENT  FOR  VARIOUS 


QUADRANT  IS  DISCRETIZED. 


1.04 


FIG.  IV-6  16  ELEMENT  CANTILEVER  BFAM 
SUBJECTED  TO  TIP  LOAD. 
NORMALIZED  TIP  DISPLACEMENT 
VERSUS  ASPECT  RATIO 


FIG.  IV- 7 16  ELEMENT  MODEL  OF  SIMPLY 

SUPPORTED  SQUARE  PLATE  SUBJECTED 
TO  UNIFORM  IOAD.  NORMALIZED 
CENTER  DISPLACEMENT  VERSUS 
ASPECT  RATIO 


D_ 

4 

FIG.  IV 


EXACT  SOLUTION 

(REISSNERS  THEORY) 

x ONE  - POINT  SHEAR 
INTEGRATION 


-8  CLAMPED  CIRCULAR  PLATE 

SUBJFCTED  TO  UNIFORM  LOAD  (48 
ELEMENT  MODEL).  COMPARISON  OF 
UNDER  INTEGRATED  SHEAR  ELEMENT 
WITH  EXACT  SOLUTION 


r/R 


FIG.  I V- 9 CLAMPED  CIRCULAR  PLATE  SUBJECTED  TO 
CONCENTRATED  LOAD  (48  ELEMENT 
MODEL).  COMPARISON  OF  UNDER- 
INTEGRATED SHEAR  ELEMENTS  WITH 
EXACT  SOLUTION 


FIG.  I V - 1 5 


UNDERINTEGRATED  PLATE 
BENDING  ELEMENT 
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